EFFICIENT ESTIMATION OF ONE-DIMENSIONAL DIFFUSION FIRST 
PASSAGE TIME DENSITIES VIA MONTE CARLO SIMULATION 



TOMOYUKI ICHIBA AND CONSTANTINOS KARDARAS 

Abstract. We propose a method for estimating first passage time densities of one-dimensional 
diffusions via Monte Carlo simulation. Our approach involves a representation of the first passage 
time density as expectation of a functional of the three-dimensional Brownian bridge. As the latter 
process can be simulated exactly, our method leads to almost unbiased estimators. Furthermore, 
since the density is estimated directly, a convergence of order 1/vTV, where N is the sample size, 
is achieved, the last being in sharp contrast to the slower non-parametric rates achieved by kernel 
smoothing of cumulative distribution functions. 

0. Introduction 

The problem of computing the distribution of the first time that a diffusion crosses a certain 
level naturally arises in many different contexts. As probably the most prevalent we mention 
quantitative finance, where first passage times are used in credit risk (times of default) as well 
as in defining exotic contingent claims (so-called barrier options). In this paper, we focus on the 
numerical computation of the probability density function of first-passage times associated with a 
general one-dimensional diffusion. 

Densities of first passage times have known analytic expressions only in very particular cases. 
The primary example is Brownian motion with certain (constant) drift and diffusion rates, where 
one uses a combination of Girsanov's theorem and the special case of standard (driftless) Brownian 
motion. The first passage time density for the latter case can be obtained by the reflexion principle 

- see, for example, [111 2.6A]. Another example of explicit form of the first-passage time density 
is the case of the radial Ornstein-Uhlenbeck process — see [7] and |10j . 

In absence of general analytic expressions for first passage time distributions, computational 
methods are indispensable and, in fact, widely used. One branch uses Volterra integral equations 

- we mention [5], [16] and [2] as representative papers dealing with this approach. Alternatively, 
one can use Monte-Carlo simulation to attack this problem. The simplest scheme uses the so-called 
Euler scheme in order to approximate the solution of the stochastic differential equation governing 
the diffusion at predetermined grid time-points ih, i = 0, 1, . . ., where h is the step size, stops at 
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the first time that the diffusion crosses the level of interest, and continuing this way obtain an 
estimator for the cumulative distribution function of the first-passage time. As the Euler scheme is 
only approximate^] it causes bias in the estimation of the probability distribution function. (See [8] 
for general discussion, [9] for the evaluation of error via partial differential equations and j3| for a 
sharp large deviation principle approach.) The issue with the bias becomes immensely more severe 
in the numerical computation of the density, as it will involve some kind of numerical differentiation 
of the (non-smooth) empirical distribution function. Even if one uses an exact simulation approach 
for the diffusion in question (which is, of course, available only in special cases), the estimator for 
the density will have huge variance. To top it all, even if all the aforementioned problems can 
be eliminated, one can never hope for convergence of the estimators to the true density in order 
1 / vN, where N is the "path-sample" size, as the problem is non-parametric. 

In this work, we offer an alternative approach which has clear advantages. First, we arrive 
at a representation of the density function in terms of expectation of a functional of a three- 
dimensional Brownian bridge. This makes it possible to estimate directly the first passage time 
density without having to rely on estimators of the cumulative distribution function, achieving 
this way the "parametric" rate of convergence l/y/N, where N is the sample size. Furthermore, 
only the three-dimensional Brownian bridge is involved in the simulation, which can be carried out 
exactly. There is an integral involving the previous three-dimensional Brownian bridge, which can 
be approximated via a Riemann sum; therefore, the error of the approximation can be estimated 
efficiently. By construction, our method significantly improves both the bias and variance of the 
density estimation obtained via the empirical distribution function. The only potential problem 
of our approach is large-time density estimation, since the thin grid that has to be used in the 
simulation of the Brownian bridge will result in high computational effort. To circumvent this 
issue, we notice that the tails of the first-passage distribution usually decrease exponentially with 
a rate that can be expressed as the principal eigenvalue of a certain Dirichlet boundary problem 
involving a second-order ordinary differential equation. This implies that a mixture of Monte-Carlo 
and ordinary differential equation techniques can be efficiently utilized and improve the quality of 
our estimator. 

The structure of the paper is simple. In Section 1, the problem is formulated and the key 
representation formula is obtained. In Section 2, we discuss the Monte-Carlo estimator of the first 
passage time density function, and study its large-sample properties. In Section 3, the relation 
between the exponential tail decay of the probability density and the eigenvalues of a Dirichlet 

■*"To account for the fact the passage can potentially happen in-between the sampled points, a Brownian bridge 
interpolation may be used — this means that the conditional probability distribution of first-passage time occurred 
in the interval given the simulated values at the end points is approximated by that of a Brownian bridge. This 
could potentially reduce the bias, but the whole scheme is still only approximate and some bias remains. 
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boundary problem is discussed. The proofs of all the results are deferred to Appendix [A] in order 
to keep the presentation smooth in the main body of the paper. 

1. A Representation of First-Passage-Time Densities 
1.1. The set-up. Consider a one-dimensional diffusion X with dynamics 

(1.1) dX t = a(X t )dt + dW t , t G K+. 

Above, IF is a standard one-dimensional Brownian motion. The restrictions on the drift function 
a that we shall impose later (Assumptions 1.2) ensure that ( |1.1[ ) has a weak solution, unique in 



the sense of probability law, for any initial condition x £ (0, oo). Let P" will denote the law on the 
canonical path-space C(M+;M) of continuous functions from M + to M that makes the coordinate 



processes behave according to (1.2) and is such that P^f-Xo = x] = 1. 

Define tq := inf {t £ R+ | Xt = 0} to be the first passage time of X at level zero. We shall 
consider the problem of finding convenient, in terms of numerical approximation using the Monte- 
Carlo simulation technique, representations of the quantity 

p a x (t) := |V« [r < t] ; x 6 (0, oo) , t € R+, 

i.e., the density of the first-passage time of the diffusion at level zero. 



Remark 1.1. The fact that we are using unit diffusion coefficient in (1.1) by no means entails loss of 



generality in our discussion. Indeed, consider a general one-dimensional diffusion Y with dynamics 
(1.2) dY t = b(Y t )dt + a(Y t )dW u t £ R+, 



where IF is a standard one-dimensional Brownian motion, such that Yq = y G M. If (1.2) has a 



weak solution unique in the sense of probability law, we may assume without loss of generality 



that a > 0. (Indeed, otherwise we replace a by |<r| in (1.2) and we obtain the same law for the 
process Y.) Consider a level £ < y. Under the mild assumption that 1/a is locally integrable, the 
transformation X = jj (1/ a{z)) dz defines a diffusion with dynamics dXt = a(Xt)dt + dWt, for 
a function a that is easily computable from b and a. With x := jf (l/a(z))dz, the first passage 
time of Y with Yq = y at level I is equal to the first passage time of X with Xq = x at level zero. 



1.2. The representation. The following assumption on the drift function in (1.1) will allow us 
to arrive at a very convenient representation for the density function p x . 

Assumption 1.2. The function a restricted on [0, oo) is continuously differentiable, and satisfies 

poo / pw \ 

/ exp I —2 / a(z)dz J dw = oo. 
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In particular, under Assumption 1.2, a is locally square integrable on [0, oo) and the function 

a 2 + a' 



(1.3) 



7 



is continuous and locally integrable. The theory of one-dimensional diffusions ensures that for 
all x G (0, oo) there exists a probability P" on C(M+,M) such that the coordinate process X has 



dynamics given by (1.1). Assumption 1.2 also ensures that X T °, which is X stopped at level zero, 



does not explode to infinity — see, for example, [TTJ, Proposition 5.32 (hi)]. 



Proposition 1.3. Suppose that Assumption 1.2 is in force. On C([0, 1]; K ), consider the prob- 



ability P BB 3 under which the coordinate process /3 is a standard 3- dimensional Brownian bridge. 
Then, 



a(v)dv I E BB 3 



exp 



uxe\ 



+ Vt/3 U 



du 



t I 7 

o 

1,0, 0) T and q x is the density given for all t G M + by 

/ i s 

- y» / • y * *J 



PS(*) = 9x(i)exp 
/joWs /or aZZ t G M+, where ei 

d.4) *W = iCW 

corresponding to the first-passage time to zero of a standard Brownian motion starting from x. 

2. Monte-Carlo Density Estimation 
We now discuss issues related to estimation of the density . For the purposes of this and the 



next section, we fix a drift function a satisfying Assumption 1.2 and we write p x for p a x in order to 
simplify notation. 

2.1. Convergence. It is clear how to get an estimate of the density p x (t) for a given t G M+, at 
least in theory. One simulates N independent paths of 3-dimensional Brownian bridge . . . , (3 N , 
and then defines the estimator p^ (t) for p x {t) via 

(2.1) £f(i) := q x (t)exp(- a(v)dv^ l^exp^-t^ 7 (juxei + Vt^|) du) , 



where recall that ^ is given in (1.4). By the strong law of large numbers, the estimator p% (t) 
converges almost surely to the true density p x (t) as N goes to infinity for each fixed t G M+. 
Moreover, the estimator p% (t) is unbiased^ and the variance of the estimator p% (t) decreases in 
the order of 1/A, for every fixed t G M+, as a direct consequence of (2.1). In order to get weak 



convergence of the whole empirical densities (p^ (i)) t£R j as well as the uniform rate of convergence 
over compact time-intervals, we introduce an additional assumption. 

2 Of course, E BB 3 [p^f (t)] — p x (t) only holds if we assume that we actually have the whole path of each Brownian 
bridge simulated exactly, which is not possible in practice. However, we can simulate exactly discretized paths of the 
Brownian bridge, and then can easily estimate the order of bias from the Riemann approximation of the integral. In 



this respect, see also [ 3.2 
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Assumption 2.1. Together with Assumption 1.2 we ask that the function 7 of (1.3) is such that 



inf 2e K + 7(2;) > —00, as well as that a' is locally Lipschitz continuous on M + with Lipschitz constant 
growing at most polynomial rate, that is, there exist constants c\ > 0, C2 > 0, n G N, such that 



(2.2) 



sup 

0<Dl <V2<K 



a'(v2) - a'jvi) 

V 2 - Vi 



< (ci + c 2 K n ) holds for all k > 0. 



The next two results are concerned with a central limit theorem for the whole density function 
estimator as well as the uniform rate of convergence on compact intervals of K4. . 



Proposition 2.2. Suppose that Assumption 2.1 holds. For NgN, define r] N := yN (p^ — p x ) ■ 

Then, the family of stochastic processes {rj N I N G N} is tight. As N — > 00 , r] N converges 
weakly to a centered Gaussian process with continuous covariance function T, where, with I(t) = 
Jq 7 (juxei + \/t/3u|) du for t G 



T(s,t) = p x (s)p x (t) exp ^ - 2 J a(w)dw^cov BB 3[exp(-s/(s)),exp(-t/(t))], (s,t) £ 



Proposition 2.3. Under Assumption 2.1, for any fixed T G R+ the sequence 



N max \f%(t) -p x (t)\ 



is bounded in probability. 



Remark 2.4. In a similar manner, for fixed x\ and x 2 in (0,oo) with x\ < X2, we may show that 



N max \pt(t)- Px (t)\ 

is bounded in probability. Moreover, under some additional conditions on differentiability of a, we 
may estimate the partial derivatives of p x (t) with respect to (x,t) via differentiating the estimator 
with respect to the variable of interest. 



3. The Rate Function 

Recall that we are dropping the qualifying "a" from "p^" in order to simplify notation. Define 
implicitly the function X x via 



Px(t) = q x (t) exp 



a(v)dv ) exp (-t\ x (t)) , t G M+. 



In other words, and in view of Proposition 1.3, we have 



(3.1) 



1 



X x (t) := —— log I Egg: 



exp ( —t 1^ 7 



uxei 



du 



for t G 
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3.1. Theoretical results. Proposition 2.3 ensures the uniform convergence of the estimator on 
finite intervals [0, T] for fixed T £ (0, oo). Of course, it is almost always the case that lim^oo p x (t) = 
0. For large t £ M + , X x (t) gives a better understanding of the behavior of the density function, as 
it represents in a certain sense the exponential decrease of p x (t); therefore, it makes more sense to 
focus on X x rather than p x . In fact, the following result implies that the function X x is frequently 
bounded on M + — this is the case, for example, when 7 is bounded from below. 



Proposition 3.1. Let Assumption 1.2 be valid. Then, 

r 7T 2 

(3.2) inf 7(2;) < inf X x (t) < limsupA^i) < inf < m(K + x) — 

where m{w) = maxo< z <„7(z) for w > 0. Furthermore, if '7 is bounded from below, then 

1 f x f 1 

(3.3) limA x (t) = — / j(u)du = / ^y(ux)du. 

40 x Jo Jo 



Remark 3.2. The inequalities (3.2) only imply bounds for the inferior and superior limit of X x (t) 
as t — > 00. In fact, it is expected that lim^oo X x (t) exists, possibly except in pathological cases. 
Let us argue below for this point on a rather loose and intuitive level. 

The stopping time tq can be approximated by the sequence (r n ) ne N, where, for all n £ N, 
7~o := inf{i > | Xt £" (0,n)} is the first exit time of X from the interval (0,n). Therefore, the 
rate function X x may be approximated by the corresponding rate functions corresponding to Tq, 
n £ N. It follows from [6] that the density function j>™ of Tq has the eigenvalue expansion (see also 
|12| . |17j for similar problems) 

8 00 
P n x it) = Wt K K < t] = ^e"^' ^(x), for (t,x) £ R + x (0,n), n £ N, 

fc=i 

for some functions {^^ | k £ N} computed from the eigenfunctions {tp^ \ k £ N} and the corre- 
sponding eigenvalues < /ij < ^ < . . . of the Dirichlet problem 

1 

(3.4) 2^"^ + a ( z )^'( z ) = ~f J,( f( z ) f° r z G (0, 7T-), 

for ip £ C 2 ' a ([0, n],R + ) with \\m. z ^Q (p(z) = = lim 2 _ >n ip(z). Thus, the limit as t — > 00 for the 
rate function of Tq is exactly the principal eigenvalue 

lim A"(i) = - lim - log ( ^l) = for all n £ N, 

<->oo t— >oo t \q. x {t)' 

which does not depend on the initial value x > 0. Since tq = linin^oo To,n, it is conjectured that 
the limit of X x (t) as t — > 00 actually exists and is equal to lim n ^.oo A thorough study of finding 
a reasonable sufficient conditions for 

lim X x (t) = lim /i™ = lim lim A"(i) 

t— >oo n— >oo n— >oo t—>oo 

to hold for x > lies beyond the scope of this paper. 
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Figure 1. Comparisons of the simulated density and rate functions in (2.1) and 
(3.6) with theoretical values for the OU process with a(z) = —z {z £ M+) and x = 1. 



3.2. Practical issues. In view of Proposition 2.2 and Proposition 2.3, the estimator (p^ (t))tgR + 



of (2.1) convergences uniformly with rate 1/yN over compact time-intervals. In practice, the 



computation of (2.1) is implemented by generating a standard 3-dimensional Brownian bridge, 



which is simulated in an exact way over a thin enough grid. The approximation error for the 



Riemann integral over the finite interval [0,1] in (2.1) can be controlled very efficiently. More 
precisely, the numerical computation of the exponential functional of the Brownian bridge in (2.1 ) 
can be carried out using the fourth-order Runge-Kutta scheme which is proposed and analyzed in 
|14| . Under appropriate mild regularity conditions on the function 7, it is shown that this numerical 
scheme is weak order of convergence 4. For this numerical issue, consult the original paper p3] 
and the related monographs [8], [13] and (15] , 

A potential problem with our estimator (p% (t))teR+ can arise for large t, that is, the density 
function at the tail. Note that what is meant here is that the relative error of the estimator of p x (t) 
tends to be large; the absolute error tends to be extremely small, as p x (t) is very close to being 
zero for large t G M+. To visualize the issue, it is helpful to study by experiment the large-time 



behavior of the rate function (3.1) when X is an Ornstein-Uhlenbeck (OU) process starting with 
x = 1. Here a(z) = —z and 7(2) = (z 2 — l)/2 for z £ K+. The first-passage time density of this 
OU process is known analytically (see, for example, PJjl equation (8)]) and reads 

1 1 n + t- coth(i) ' 



(3.5) 



Pi(t) 



2tt sinh 3 / 2 (t) 



exp 



for t G R+. 
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The true density (|3.5[) and the estimated density (2.1 ) with N = 100 simulations are shown over 



interval [0, 10] in Figure [ija) . Note that even with this small number of simulations, the two curves 
are almost indistinguishable. Figure [T](b) contains graphs of the more refined rate functions. In 
this scale, one can see that the estimation of the exponential rate of decay of the density for large 
t is not as good. (As can be seen from[I](b), on the interval [0,7], the true and estimated rate 
functions almost coincide. However, on the interval [7,20], the estimated rate functions are 



much larger than the true rate function X x , which implies that the estimator (2.1) underestimates 



the tail probability.) In fact, the estimated asymptotic rate seems to increase linearly instead of 
converging to a finite limit. This becomes clear once one notes that, in this OU example, 



-t J y(uxe l + Vt(3 u ^Jdu = t 2 ^J^ \p u \ 2 du + t 3/2 x 



(ei,p u ) du + 1 : 
io 6 



therefore, the estimator for A^ for sample size A'gN becomes 



(3.6) 



x 
~6~ 



x I (ei,0l) du 



Observe that the leading term in the estimator of X x will be increasing linearly in t. 

In order to overcome this poor situation in the tail, one can use a mixture method combining 



the estimator (2.1 ) on the finite interval [0, T] and an estimator for the tail probability of the form 

c*q x (t) exp(— Xt), for t S [T, oo), 
for some choice of large threshold T, where A is the principal eigenvalue of the Dirichlet prob- 



lem (3.4) for some choice of large threshold n, and c* is chosen so that the density estimator 



is continuous. The principal eigenvalue can be numerically computed from the Sturm-Liouville 
problem 



exp { 2 I a(u)du) ip'{z) 



2A< x|> ( 2 I a(u)du)(p(z) for z € [0, n] 



(see [20 1) with the same Dirichlet boundary conditions, by use of either the variational method or 
the Liouville transform method. In the variational method, the principal eigenvalue is obtained 
by minimizing numerically the corresponding Rayleigh quotient. The Liouville transform method 
turns the Sturm-Liouville equation into a Schrodinger-type equation, which is then numerically 
solvable with discrete approximation. Both numerical methods are well studied — see [1] and their 
references within. 



Appendix A. Proofs 



A.l. Proof of Proposition 1.3, Consider the non-negative P"~ su P ermar t m g a l e 

1 



(A.l) 



exp 



At 



At 

- V .Ml - 2 J^ a 2 (X u )du 
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As follows from a modification of [HJ Exercise 5.5.38] for the restricted state space [0, oo), Z is a 
P"-martingale. With <Q X being the probability on C(M+;M) that makes the coordinate process X 
behave like a Brownian motion staring from x and stopped when it reaches level zero, Girsanov's 
theorem implies that 



(A.2) 



dP£ 



Zt. 



for all t E K+. Moreover, since X T0 = 0, Ito's formula (under Q x ) implies that, on the set {tq < oo}, 



/ a{v)dv = / a(v)dv = / a(X u )dX u + - / a\X t 

JO Jx JO * Jo 



)du. 



Combining this with (1.1) and (1.3), the stochastic exponential defined in (A.l) under Q x satisfies 

fTO 



1 



■"TO 



exp / a{X u )[dX u - a{X u )du] + - / a~{X u )du 



o 



o 



(A.3) 



exp 



exp 



exp 



TO 



a(X u 



a(v)dv 



1 f T0 



+ 



i n , 



a (X u )du 
1 



a (X u )du 



TO 



a 2 (X u )du 



f a(v)dv - [ j(X u )du), 
o Jo ' 



on {tq < oo}, where W^ x is a standard Brownian motion under <Q X . Therefore, (A.2) and (A.3) 
imply 



^[r <t]=E c 



1 



Z 



*At 



L {T0<t} 



exp 



o 



fTO . 

a{v)dv-j j(X u )du)I 



{T0<t} 



for (t,x) E (M+) 2 . 

Note that the density function of tq under Q x is given by q x in (1.4). Using the regular conditional 



-expectation of exp(J Q T0 j(X u )du), given tq = t, we can write 



(A.4) p a x (t)=F a x [T £dt]=q x (t)exp 



a(v)dv)E Qx 



exp 



j(X u )du 



t = t 



Given tq = t, the regular conditional (Ql^-distribution of (Xt_ s , < s < t) is that of a 3-dimensional 
Bessel bridge from to x over [0, t] as a consequence of |19} Proposition VI. 3. 10 and Proposi- 
tion VII. 4. 8]. On the canonical space (C([0, 1],M 3 ),P BB 3) with coordinate process /3, the process 
{\(s/t)xei + y/i/3 s u\, < s < t} has the exact law of the aforementioned Bessel bridge. Therefore, 



exp 



TO 



j(X u )du 



T = t 



E BB 3 



exp 



-xei 



+ Vtp, 



s/t 



ds 



exp -i / 7 



uxei 



du 



; (t,x)e(R + f 



Combining this with (A.4), the proof of Proposition 1.3 is complete. 
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A. 2. Proof of Proposition 2.2 , The following technical result is the backbone of the proof. 



Lemma A.l. Suppose that Assumption \2.1\ holds, and define 
(A.5) £(t) := exp(-t/(t)) := exp ( - t J 7 



uxei 



+ Vi(3 u 



du), t £ 



Then, we have 

!£(*)- < ®T\t-s\ for alls G [0,T], and t G [0,T], 
w/iere E BB 3 [|$r| m ] < oo for all Tel + and m G N. 

Proof. First, note that (2.2) in Assumption |2 . 1| implies that for every k > 0, 

< k'(0)| + |o'(ui) - a' (0)| < | a' (0)| + ci«i + C2/s n vi; < «i < «, 
|a(«)|= a(0)+ / a'(u)dit < |a(0)| + / |a' (u)|du < |o(0)| + |o'(0)|« + cik 2 + c 2 K n+2 , 



|a(«i) - a(y 2 )\ < / o'(u)du < (|a'(0)| + cik + c 2 K n+1 )\v 2 - vi\; < v t < v 2 < k. 

J Vl 

Using these inequalities, we obtain estimates for 7 for every k > and every < v\ < v 2 < k, 



h(vi)\ 



a 2 {vi) + a'(y 1 ) 



1 



< ^(|a(0)|+|a / (0)| K +ci K 2 + C2K n+2 ) 2 +^(|a , (0)|+ C i K +c 2K n+1 ) =: ^(k), 



h{vi) - j(v 2 )\ < \ \a 2 {vi) - a 2 (v 2 )\ + - \a'(yi) - a'(v 2 )\ 



< - \a(vi) + a(v 2 )\ \a(vi) - a(v 2 )\ + -(ci + c 2 K n ) |t>i - v 2 \ 
(A.6) 2 2 

< "(|a(0)| + |a'(0)|« + cik 2 + c 2 K n+2 ) (|a'(0)| + c x k + c 2 K n+l ) + i( Cl + c 2 K n )] \ Vl - v 2 \ 

=: (p 2 (n) \vi - v 2 \ , 

where M + 3kh> J = 1) 2 > are polynomial functions of k G M+ and do not depend on v\, v 2 . 

Fix T G M+. For s G [0, T], t G [0, T] and u G [0,1], consider the random variables K = 
\/Tmaxo< n <i|/3 n | + x, v(s,u) = \uxe\ + y/s/3 u \, v(t,u) = \uxe± + Vi(3 u \. Using the estimates 
established before, we obtain estimates for I(t) in (A.5): 

\I(t)\= [ -f(\ux ei + Vi(3 u \)du < [ \-f(\ux ei + Vif3 u \)\du = [ \j(v(t, u))\du < </>i(k), 
Jo Jo Jo 

s ll 1 (t) — I 1 (s)\ < s / \j(v(t,u)) — j(v(s,u))\du < s / \v(t,u) — v(s,u)\(p 2 (n)du 



(A.7) 



<8(Vt-Vs) f \/3 u \<p 2 (K)du < S j* s \ ^-=,y 2 (n) 
Jo Vt + vs v T 



2VsT 
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where we have used (A. 6) in the second inequality, since < v(s,u) < k, < v(t,u) < k, used 
\v(t, u) —v(s,u)\ < {y/t— \fs)\(3 u \ in the third inequality, and maxo<«<i|/3 M | < k/VT in the fourth 
inequality for < s < t < T. 

Finally, since 7 is bounded from below by Assumption 



2.1 



so is / in AAM, that is, I 1 ^) > 
inf 2e iR + 7(2;) > —00 for every t > 0. With this observation, because of monotonicity and differen- 
tiability of exponential function, we obtain for < s < t < T, 



-U(t) 



-al{a) 



< ( e -Tinf 7 v ^ | t/ ^ _ s/(s) 



e w — e 

= c 3 \(t - s)I(t) + s(I(t) - I(s))\ <c 3 \t- s\ \I(t)\ + c 3 s \I(t) - I(s)\ , 

where C3 := exp(— Tinf 7) V 1 < 00. Combining this with the estimates for |/(t)| and s\I(t) — I(s)\ 
in (A. 7), we obtain, for < s < t < T, 



\m-t(s)\ <c 3 (W«)+ 



where tp\ and ip 2 are defined in (A.6), and hence (^3 can be written as a polynomial function of k 
whose coefficients do not depend on s nor on t but on T. Letting be 993 (k), and noting that 
all positive integer moments of maximum of standard 3-dimensional Bessel Bridge are finite ([181 
Corollary 7]), we conclude the proof of Lemma A.l □ 



Let us define u{t) := E BB s for < t < T. It follows from Lemma A.l that £ is locally 
Lipschitz continuous and moreover, 



E 



BB 



ie(*)-*oo-("(t)-"W)i ! 



< C4 \t — s\ 



< s < t < T. 



Since the random paths = 1,...,N} are independent and identically distributed, for any 

s, t < T we obtain 



E 



ri N (t) 



BB 



V N (s) 



(A. 



< 



(q x {t)) exp - 2 
max q x (t)) 2 exp ( — 2 



a(«)dw)E BB 3 [|e(t) - ^(s) - (v(t) - v{s))\ 



a(v)dv )c4 |t 



- c 5 |t 



v 0<t<T 

where the constant C5 depends on T, x but not on N,s,t. This inequality is a sufficient condition 
for the tightness of the sequence {r] N \ JVgN} of continuous stochastic processes starting at 
in C(IR+,]R) — see [111 Problem 2.4.12]. By the usual multi-dimensional central limit theorem, 
for each n > 1, < h < ■ ■ ■ < t n < 00 the sequence { (^(ii), 77^*2), • • -,ri N (t n )) \ N £ N} of 
random vectors converges in distribution to a Gaussian random vector with mean zero and and 
variance-covariance matrix (r(ij, tj)) 1<{ - <n , where 

T(s,t) = exp ( - 2 y a(v)d^cov BB 3[exp(-s/(s)),exp(-tJ(t))], (s,t) G 
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Therefore, we conclude that the tight sequence {rj N | N G N} converges weakly to a continuous 
Gaussian process with mean zero and continuous covariance function V. 



A. 3. Proof of Proposition 2.3 Define that Gaussian tail function via 

e" 



dy, for zel. 



Furthermore, for fixed T G M+, define the modulus of continuity in L 2 : 

\2lV2 



(A.9) 



Mh) := , w max [E BB 3 (»/(*) - r/( S )) 2 ] 7 , for /i G [0,T]. 



It follows from (A. 8) that V"r(^) < for h G M+; therefore, ^(e )du < oo. We recall 



the following Fernique's inequality for Gaussian processes, which we shall use. 



Lemma A. 2 (Fernique's inequality — see (2.2) of [1]). If the function ipj- in (A.9) satisfies 
V'r(e _ ^ 2 )dy ; then for that fixed T > and any integer m > 2, 



max \n(t)\ > Ci(T,m)z 

0<t<T 



< C 2 (m)$(z), for all z > (1 + 41ogm) 1/2 , 

where Ci(T,m) := max (Sit)e[0 T]2 T(s,t) 1 / 2 + (2 + ^2) ip T (Tm~ y2 )dy, and C 2 := 5m 2 v / 2^/2. 
The weak convergence of {ry^ | A G N} to r\ and the invariance principle for the maximum 



function imply that the sequence of normalized maxima f y Amax tg [ ,T] |Pir (*) ~~ Paj(*)| ) con ~ 
verges weakly to maxo<t<T ^(^), where r) is the limiting Gaussian process, as N goes to infin- 
ity. Since the law of the last random variable does not charge oo, we conclude that the family 
j\/iV max^^T] \f)x it) — p x {t)\ | N G n| is bounded in probability. 

A. 4. Proof of Proposition |3?L For the lower bound in (3.2) let us observe 



logE BB 3 



exp \-t / 7 



by Assumption 2.1 Therefore, 

1 



inf 7(2) < inf 



f 



logE BB a 



uxei 



exp -t / 7 



du 



< -t inf 7(2;); i G 



uxei 



du 



inf X x (t). 



For the upper bound in (3.2), for a fixed k > let us consider A K := |vi maxo< u <i|/3 u | < k\. On 
^4 K , \uxe± + \fij3 u \ < k + x holds for < u < 1; hence, 

Jq V / 0<2<K+a: 

where m(w) := maxo< z < 1 _„ 7(2;) for u; > 0. It follows that, for t G K+, 



E 



BB 



exp -i / 7 



uxei 



du 



> E 



BB 



exp -i / 7 



uxei 



du Ia k 



> exp (—tm(K + x)) P BB 3 yr max |/3 U | < k 

L o<m<i 
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and hence 
(A.10) 



1 



X x (t) < m(n + x) log 



t 



BB a 



max \B U \ < nt 

0<u<l 



-1/2 



The distribution for the maximum of the absolute value of the standard 3-dimensional Bessel bridge 
|/3 1 is known — see, for example, |181 (5)]. More precisely, we have 

2 [2t A nvr 



(A.11) 



BB' 5 



max \B U \ < Kt 
0<u<l 



-1/2 



* J 3/2( 



exp 



7r 2 n 2 



t ) , for k > and i 6 



where J3/2 is the Bessel function of index 3/2. In particular 

1 



lim 

t— K>0 



t 



BB 



3 log 



max \B U \ < nt~ 

0<«<1 



1/2 



7T 

2^' 



as only the first term in the summand in the series in ( A. 11 ) will play a role in the limit. Combining 



the last limiting relationship with inequality (A. 10), we obtain 

lim sup X x (t) < m\K + x) + 



Upon minimizing the right-hand side of the last inequality, the upper bound in (3.2) is obtained. 



Finally, to verify (3.3), observe that 7 being bounded from below implies that the random vari- 
ables exp(— t Jq j(\uxei + y/i/3 u \)du) are uniformly bounded for small t G M. + . Then, de L'Hospital's 
rule and the bounded convergence theorem give 



limA x 0) = (-l)lim 



40 



40 



(d/dt)E BB 3[exp(-t J 7(|iiacei + VtB u \)du)] 
E BB3 [exp(-t 7(|nxei + Vi/3 u \)du)] 



'j(ux)du 



1 

x Jo 



j(u)du. 
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